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Abstract 
(N 

Relativistic Continuum Random Phase Approximation (CRPA) is 
c| used to investigate collective excitation phenomena in several spher- 

ical nuclei along the periodic table. We start from relativistic mean 
q field calculations based on a covariant density functional with density 

dependent zero range forces. From the same functional an effective 

1 ^ i interaction is obtained as the second derivative with respect to the 

density. This interaction is used in relativistic continuum-RPA cal- 
culations for the investigation of isoscalar monopole, isovector dipole 

CnJ and isoscalar quadrupole resonances of spherical nuclei. In particu- 

lar we study the low-lying El strength in the vicinity of the neutron 
evaporation threshold. The properties of the resonances, such as cen- 
troid energies and strengths distributions are compared with results 

<^ of discrete RPA calculations for the same model as well as with ex- 

perimental data. 

PACS numbers: 21.30.Fe, 21.60. Jz, 21.65.+f, 21.10.-k 
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1 Introduction 

Density functional theory (DFT) provides a very successful description of 
nuclei all over the periodic table. Based on relatively simple functionals, 
which are adjusted in a phenomenological way to the properties of infinite 
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nuclear matter and a few finite nuclei, this theory allows a highly accurate 
reproduction of many nuclear structure data, such as binding energies, radii, 
deformation parameters of finite nuclei and their dependence on mass number 
and isospin. In addition to these static properties, one can use the nuclear 
response to external multipole fields to investigate the dynamics of such 
systems. In the framework of time dependent density functional theory, this 
response can be calculated from the linearized Bethe Salpeter equation using 
an effective interaction derived from the same functional. 

A very successful scheme of this type is covariant density functional theory 
(CDFT). It is based on Lorentz invariance, connecting in a consistent way the 
spin and spatial degrees of freedom of the nucleus. Therefore, it needs only 
a relatively small number of parameters which are adjusted to reproduce 
a set of bulk properties of spherical closed-shell nuclei. Numerous works 
have shown that observations involving both ground state and excited state 
phenomena, can be nicely interpreted in a relativistic framework. 

The most popular applications of this type are based on the Walecka 
model P, where the nucleus is described as a system of Dirac nucleons inter- 
acting with each other via the exchange of virtual mesons with finite mass and 
the electromagnetic field through an effective Lagrangian. In the mean field 
approximation this yields to various contributions to the nuclear self energy 
depending on the quantum numbers of these mesons. Early investigations 
have shown that this simple ansatz is not able to describe the incompress- 
ibility of infinite nuclear matter nor the surface properties of finite nuclei 
such as nuclear deformations. For that reason, a medium dependence has 
been introduced by including nonlinear meson self-interaction terms in the 
Lagrangian [2J. 

Several very successful phenomenological RMF interactions of this type 
have been adopted, as for instance the popular set NL3 [3]. Closer to the 
concept of density functional theory are models with an explicit density de- 
pendence for the meson nucleon couplings. This density dependence can be 
calculated from first principles in a microscopic Dirac-Brueckner scheme jl] 
or it can be adjusted in a completely phenomenological way to properties of 
finite nuclei [HI E] . 

One of the advantages of density functional theory is the fact that with a 
proper choice of the parameters the success of RMF for nuclear ground states 
ensures also a good basis upon which one can apply time-dependent density 
functional theory to study nuclear excitations. In order to investigate the 
dynamic behavior of the nuclear system, one considers oscillations around 
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the self-consistent static solution. This can be done by solving the time 
dependent relativistic mean field equations (TDRMF) [7] or, in the limit 
of small amplitudes, by using the relativistic random phase approximation 
(RRPA) [8]. The corresponding eigen modes can be determined either by 
diagonalizing the RRPA equation in an appropriate basis or by solving the 
linear response equations in a time-dependent external field. This requires a 
matrix inversion for given frequency u. 

These two methods lead in principle to exactly identical results. There 
are, however, cases where one of them is clearly preferable. The proper 
treatment of the coupling to the continuum is such which can be 

solved in a very elegant way, by the solution of the Bethe Salpeter equation 
within the response formalism. 

We recall that the spectrum of the Dirac equations has a discrete and a 
continuous part. For the ground state properties of the nucleus, one needs 
only the single particle wave functions of the occupied orbitals in the Fermi 
sea. They are either determined by solving the corresponding differential 
equations in r-space or by expansion in an appropriate basis, given for for 
instance by a finite number of eigenfunctions of a harmonic oscillator [9] 
or of a Saxon- Woods potential in a finite box [ID] . For the bound states 
both methods yield the same solutions with high accuracy. However, this is 
no longer true for the states in the continuum. Here we have, in the first 
case scattering solutions in r-space for each energy with proper boundary 
conditions, while in the second finite number of discrete eigenstates 

which depend strongly on the dimension of the expansion. They provide only 
a basis and have little to do with physics. 

These discrete eigenstates lead to a finite number of p/z-configurations for 
the solution of response equations, with a discrete spectrum. They provide us 
with the so called spectral representation of the response function in contrast 
to the continuum representation, where the exact scattering states with the 
proper boundary conditions are used at each energy. 

Self-consistent relativistic RPA (RRPA) calculations have a long history. 
The early investigations in the eighties [ITJ [131 [HJ EES] were based on the 
Walecka model with linear meson-nucleons couplings. They were able to 
describe the low-lying negative-parity excitations in 16 by the method of 
matrix diagonalization [11] , isoscalar giant resonances in light and medium 
nuclei [13] by the solution of the linear response equation in the spectral rep- 
resentation, and the longitudinal response for quasi-elastic electron scattering 
with a proper treatment of the continuum. 
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The first RRPA calculations based on non-linear models were carried out 
in the spectral representation including only normal particle-hole (ph) pairs 
with particles above the Fermi energy and holes in the Fermi sea. This 
seemed to be a reasonable approximation, since the configurations formed 
by particles in the Dirac sea and holes in the Fermi sea (a/i-pairs) are more 
than 1.2 GeV away from the normal p/i-pairs. Indeed, a proper coupling 
to the Dirac sea and current conservation was neglected in these investiga- 
tions. They showed considerable deviations from the results obtained form 
time-dependent RMF-calculations with the same Lagrangian, particularly for 
isoscalar excitations PS] . A fully self-consistent treatment with current con- 
servation requires the inclusion of a very large number of a/i-pairs connected 
with a considerable numerical effort. Most of the very successful applications 
of RRPA theory based on non-linear meson-nucleon coupling models in the 
last ten years have been carried out in this way [T7J HH HH [2Ql [21] . 

There are also relativistic continuum RPA calculations based on the non- 
spectral representation of the response function using the single particle 
Green's function in the continuum with proper boundary conditions [14J. 
These calculations are done for meson exchange forces with finite range. The 
early investigations were based on linear models. Later on the method was 
generalized to include non-linear coupling terms between the mesons [22J. 
This leads to a a more sophisticated density dependence which is crucial for 
a realistic description of giant resonances in nuclei [221 [23] • 

Of course, because of the finite range of the effective force these mod- 
els are relatively complicated not only for static applications to triaxially 
deformed or rotating nuclei, but also for investigations of nuclear dynam- 
ics, such as the solution of the relativistic RPA or linear response equations 
for the description of excited states. In particular one needs simpler forces 
for applications going beyond the mean field approach such as Particle Vi- 
brational Coupling (PVC) [21] or configuration mixing calculations in the 
framework of the Generator Coordinate Method (GCM) [25]. Therefore over 
the years several attempts have been made to develop relativistic point cou- 
pling (PC) models with forces of zero range [2E], in analogy to non-relativistic 
Skyrme-functionals. but only recently parameter sets have been found, which 
are comparable in quality to the density dependent meson-exchange mod- 
els [23128]. 

PC models do not contain mesonic degrees of freedom and are therefore 
closer to the philosophy of the density functional theory. Their essential ad- 
vantage is of course the fact that the zero range of the effective interaction 
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reduces considerably the numerical effort in practical applications. Because 
of their simplicity they are nowadays much used in many complex calcula- 
tions going beyond the mean field approach [2U [25] . However, so far they 
have not been used much for the dynamic investigations and it is only quite 
recently that a code has been developed to diagonalize the RPA equations for 
relativistic Point Coupling models [21] and it has been shown that this lat- 
ter approach reproduces excitation and collective phenomena, in particular 
Giant Multipole Resonances, with a quality comparable to that of standard 
finite-range forces. 

This manuscript is devoted to an investigation of relativistic point cou- 
pling models with an exact treatment of the coupling to the continuum. The 
relativistic response equations are solved both in the continuum and in the 
spectral representation and the corresponding results are compared. We use 
the Lagrangian PC-F1 [27J, which is capable of reproducing a wide range of 
experimental data. 

The paper is organized in the following way: In Sec. [2] we present the 
main characteristics of the point coupling RMF theory, while the relativistic 
RPA equations are derived in Sec. |3j The proper treatment of the continuum 
in connection with point coupling models is discussed in Sec. [4] and in Sec. [5] 
we finally present applications of this method for the spectra of in spheri- 
cal nuclei. In particular we calculate the strength function of Isoscalar and 
Isovector Giant Resonances as well as their contributions to their respective 
energy weighted sum rules. The results are summarized in Sec. [6] 

2 Relativistic mean field theory of zero range. 

As in all the relativistic models, the nucleons are described as point like 
Dirac particles. In contrast to the Walecka model, however, where these 
particles interact by the exchange of effective mesons with finite mass, point 
coupling models [26] neglect mesonic degrees of freedom and consider only 
interactions with zero range. In principle, these models are similar to the 
Nambu Jona-Lasinio model [30] used extensively in hadron physics. There is, 
however, an important difference: in order to obtain a satisfactory description 
of the nuclear surface properties one needs gradient terms in the Lagrangian 
simulating a finite range of the interaction. 

A general point-coupling effective Lagrangian is constructed to be con- 
sistent with the underlying symmetries of QCD (e.g., Lorentz covariance, 
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gauge invariance, and chiral symmetry). It should in principle contain every 
possible term, allowed by these symmetries, but at the same time should also 
be described by the least possible number of parameters in order to give a 
quantitative solution. 

In this work we use the point coupling Lagrangian introduced by Buer- 
venich et al. in Ref. [27]. It presents an expansion in powers of the nucleon 
scalar, vector and isovector-vector densities. The Lagrangian 

C = £frcc + Af + Aiot + Aer + Am (1) 

consists of the term for free nucleons: 

Aree = ^(ij^ - m N )^p, (2) 

the term for normal four-fermion interactions 

Af = - y (W) (W) - °y (h^) ( 3 ) 

the term for higher order terms leading in mean field approximation to a 
density dependence 

Act = -y (« 3 - f (W 4 - ^[(^)(^)] 2 , (4) 

the term containing derivative terms which simulate in a simple way the finite 
range of the forces: 



- - ^(d^T^id^T^) (5) 

and finally the electro-magnetic part 

Am = ~V - |(1 - r 3 )^ 7 ^- (6) 

In these equations, t/> represents the nucleon spinors. The subscripts S and V 
are attributed to scalar and vector fields, while the subscript T is attributed 
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to isovector fields. As usual, vectors in isospin space are denoted by ar- 
rows, where symbols in bold indicate vectors in ordinary three-dimensional 
coordinate space. 

From this Lagrangian and the corresponding energy momentum tensor 
we can derive a relativistic energy density functional. It has the form: 

WLM = J d 3 rH(r,t), (7) 

where the energy density 

H(r, t) = H kin (r, t) + H int (r, t) + H cm (r, t) (8) 
consists of a kinetic part 

A 

#kinO, t) = ^ 4>i(r, t) (ctp + (5m -m) ^(r, t), (9) 

i 

an interaction part 

+ °f Uf + y Cm") 2 + y W (io) 

+ ~^3tV ■ JTVn + ~YJtV ■ A (jTV)n 

and an electromagnetic part 

H em (r, t) = \f, v F^ - F°»d A„ + eA,f p . (11) 
The interaction part depends on the local densities: 

A 

p s (r,t) = Y,Mr,t)Mr,t), (12) 

i 

A 

Pv(r,t) = ^2^i(r,t)-f ipi(r,t), (13) 

i 

A 

Prs(r,t) = ^2^i(r,t)T 3 ^i(r,t), (14) 

i 

A 

pTv(r,t) = ^^(r,t)r 3 7o?Mr,t) (15) 
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and currents 



j£(r,t) = ]T^M)7^(M, (16) 

i 

% v (r,t) = Y^UrJW^t)- (17) 

i 

As in all relativistic mean field models, the no-sea approximation is used in 
the calculations of the nuclear densities by summing only over the single- 
particle states with energies in the Fermi sea. Vacuum polarization effects 
are not taken into account explicitly but only in a global way by the correct 
choice of the Lagrangian parameters. All interactions in the Lagrangian (fl]) 
are then expressed in terms of the corresponding local densities 

Many effects, which go beyond mean field, seem to be neglected on the 
classical level, such as Fock-terms, vacuum polarization, short range Brueck- 
ner correlations etc. However, the coupling constants of the method are 
adjusted to experimental data, which, of course, contain all these effects and 
many more. Therefore these effects are not neglected. On the contrary, they 
are taken into account in an effective way. This concept of RMF methods is 
therefore equivalent to that of density functional theory. 

The time-dependent variational principle 

6 J {^Ij^it)) ~ E[p(t)\] dt = (18) 

allows us to derive from the energy density functional E[p] an equation of 
motion for the time-dependent relativistic single particle density: 

A 

p(r,r',t) = Y,\Mr,t))(Mr'M (19) 

i 

which has the form 

id t p(t) = [h(p(t))r P (t)}. (20) 

The self energy, i.e. the single particle hamiltonian h(p(t)) is obtained as 
the functional derivative of the energy density functional with respect to the 
relativistic density matrix: 

ft=^P. (21) 
dp 
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This yields the Dirac hamiltonian: 

h = a[-iV - V(r, t)] + V(r, t) + (3{m + S(r, t)) (22) 

with the self-consistent scalar and vector potentials 

S(r,t) = E s (r,t)+i ! -2 TS (r,t), (23) 
V i (r,t) = E»(r,t)+T-E£(r,t). (24) 

The nucleon isoscalar-scalar, isovector-scalar, isoscalar-vector and isovector- 
vector self-energies are density dependent and defined by the following rela- 
tions: 

£s = a s ps + PsPs + IsPs - SsAps, (25) 
Y,ts — ®tsPts — Sts^-Pts, (26) 
^ = a v p v + ivPv ~ <VAp y - eA» 1 -^, (27) 
Ej. = cxtvPtv — Stv^Ptv- (28) 

Here we have neglected retardation effects, i.e. second derivatives with re- 
spect to the time for the various densities. 
In the static limit we have 

[h(p),p]=0, (29) 

thus the static density po is obtained from the solution of the self-consistent 
Dirac equations upon all the nucleons with eigenvalues Ek and eigenfunctions 

h\Mr)) =£k\Mr))- (30) 
For spherical symmetry the spinors have the form: 

l *"" (r)> = ;li*.(r)y. m (n)J- (31) 

The subscripts n, k and m are principal and angular momentum quantum 
numbers; k — ^(j + |) for j — I ± |, where j and I are the total and the 
orbital angular momenta of the nucleon. As usual, m is the z component 
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Coupling const. 


PC-F1 


a s 


-14.935894 


Ss 


-0.634576 


ay 


10.098025 


Sv 


-0.180746 


a T s 


0.0 




0.0 


axv 


1.350268 


OTV 


-0.063680 


Ps 


22.994736 


Is 


-66.769116 


Iv 


-8.917323 



Table 1: The coupling constants in the parameter set PC-F1 resulting from 
the fitting procedure in Ref . [27] . The units are [fm -2 ] for the constants a of 
the quadratic terms, [fm -4 ] for the constants S of the derivative terms, [fm -5 ] 
for the constants (3 of the cubic terms, and [fm -8 ] for the constants 7 of the 
quartic terms in the Lagrangian. 



of the total angular momentum. The spherical spinors y Km (fi) are given in 
terms of spherical harmonics Y[ mi (ft) and Pauli spinors Xm s as: 

y K m(ty = ^2(^m s lmi\jm)Y lmi (^)Xms, ( 32 ) 

minis 

while the functions fi(r) and gi(r) satisfy the static radial Dirac equations: 

y + S ~d r + ^ \ffi(r)\_ff i (r)\ 
d r + ^ V-S-2m J \ 9l (r) J ~ \ 9i (r) ) ^ [66) 

The point coupling Lagrangian used in this work contains eleven cou- 
pling constants. Based on an extensive multi parameter \ 2 minimization 



10 



procedure, Biirvenich et al. [2Z] have adjusted the parameter set PC-F1 to 
reproduce ground state properties of infinite nuclear matter and spherical 
doubly closed shell nuclei. This set is listed in Table [T] and it has been tested 
in the calculation of many ground state properties of spherical and deformed 
nuclei all over the periodic table. The results are very well comparable with 
reasonable effective meson-exchange interactions. 

The nuclear ground state is defined as the equilibrium point of the func- 
tional ([7]), thus, is associated with the density which minimizes £?rmf[p]- Fur- 
thermore, small oscillations around this equilibrium point correspond to the 
vibrational nuclear states. They are usually described within the harmonic 
approximation, that is, using linear response theory. In nuclear physics, this 
is the so called Random Phase Approximation (RPA) which has been already 
mentioned in our discussion and will be described in more detail in the next 
section. 



3 Relativistic RPA formalism 

Under the influence of an external field F(uj) oscillating with the frequency 
uj the nucleus is excited. The cross section of this process is proportional to 
the strength function: 

S(uj) — Im F a p*R a p a ii3i{uj)F a ifji 

71 * — ' 

a/3a'/3' 

:= --ImR FF (u), (34) 

7T 

where F a p is the operator inducing the reaction and R a pys(w) is the response 
function which, in an arbitrary representation indicated by the Greek indices 
a, P, . .. (e.g. the (r, ^-representation) is defined as: 

d r\ ( (°H a <x\ u )( u \ a a' a P f \°) (v^a^ty^a+a^v) \ 

iWH-l.j u-Ev + Eo + iri u + Ev -E + i V /' (35) 

The imaginary part if] is infinitesimal and is introduced in order to fulfill the 
proper boundary conditions and to prevent R(u) from diverging at u = E v — 
E . We use here the response derived from the retarded Green's functions 
as defined in Ref. 
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In the independent particle model, |0) is the Slater determinant of the 
ground state, formed by the self-consistent solutions of the Dirac equa- 
tion (30) and \v) = a^a/jO) are p/z-states, while E and E u are the corre- 
sponding energies. In the basis \k), where the single particle hamiltonian (22) 
is diagonal we obtain the free response function: 



R 



n k - m 



u - e k + ei + irj 



with the occupation factors: 
n k = (0\a^a k \0) = 



1 for hole states with e k < Ef 
for particle states with e k > Ef 



The full response of Eq. (35) contains the transition densities: 

Pa/3 = (0\a+a a \u). 



(36) 



(37) 



(38) 



They can be deduced from the time-dependent density matrix in Eq. ([T9j) , 
which is derived from the variational principle in Eq. (18). 

In the small amplitude limit one uses the linear response approximation 
to obtain the full response R(u) of Eq. (35) as the solution of the linearized 
Bethe-Salpeter equation: 

R af 3 a 'p>(u}) = R°p alf3 ,(uj) + ^2 R^si^V^ys'Rys'a'p'iu). (39) 

The relativistic residual interaction is found as the second derivative of the 
energy density functional ([TJ) with respect to the density matrix 

h _ PE\p] 

Once again, we have neglected retardation and this effective interaction has 
to be calculated at the static density. 

In a short hand notation the response equation (39) has the formal solu- 

(41) 



tion 



R(u) = (1 - i? V)^ ph ) _1 # V) 
or introducing the inverse of R° we have 

1 



R(u) 



RHuj)- 1 - VP h 



(42) 
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The evaluation of the strength function (34) requires therefore three steps. 
The starting point is the calculation of the free response function R°(u>). In 
the next step one determines the interaction V ph and finally one solves the 



response equation by the inversion (41). In details there are several methods 
to proceed. In particular one can choose various basis sets to solve these 
equations. 

a) As we have seen in Eq. (36) the free response has a particularly 



simple form in the basis of Dirac spinors (Dirac basis) diagonalizing the self- 
consistent mean field equation (30). This is in particular simple for cases 



where the Dirac equation is solved in a discrete basis, as for instance the 
oscillator basis [H] or in a Saxon Woods basis [TU] determined by the solution 
of the Dirac equation in a box with finite size. However, the simplicity in the 
calculation of R°(uj) is compensated by the computational effort required in 
the next steps. First we have to calculate a large number of matrix elements 



for the interaction (40) in the basis of the corresponding p/i-states and in 
a second step the matrix (1 — R°(u)V ph ) has to be inverted for each value 
of the frequency u. In general the number of single particle states is rather 
large and this leads to a huge number of p/i-states, requiring considerable 
computational sources, not only in memory but also in computer time. This 
is in particular a problem in the case of deformed nuclei. By this reason this 
method can only be used successfully for light spherical nuclei, where the 
number of p/i-states is limited. 

b) The inversion is particular simple in the RPA-basis. Inserting expres- 



sion (36) into Eq. (42) we find that the response function is equivalent to the 



resolvent of the RPA matrix 



where 



R°(u) 



Aphp'h' 

Bphp'h' 



A B 

-B* -A* 



(e p - 6 h )5p p/ 8 hh i + V^, h , 



T/Ph 
v phh'p' 



(43) 

(44) 
(45) 



Of course, the calculation of this matrix requires the same numerical effort 
as the evaluation of V ph in the Dirac basis discussed above. However there 
exist standard routines for the diagonalization of the RPA-matrix 



A 

-B* 



B 

-A* 



X 
Y 



X 
Y 



(46) 
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and this diagonalization has to be carried out only once, whereas the inversion 
of the response equation has to be done for each value of the frequency u. In 
the RPA-basis given by the eigenvectors \fi) the reduced response function 
defined in Eq. (52) has a particular simple form 

\Qt\ri(y\Qc>\0) (p 

lu - fi M + it] uj + fi M + it] 



cc'{ ) ,,, _0.. 4-7'n i.i 4- O 4- in ' ^ ' 



H>0 



Using 



(0\F\ri=Y,F»(XZ> + Yik) (48) 

ph 



we find for R F p(uj) 



and for the strength function in Eq. (34) 



5(o; + z-) = --lmR FF (u + i—) (50) 

2 7T 2 



A 



27T (^-O )2 + 1 A 2 



Here A is a smearing parameter, which introduces a folding with a Lorentzian 
and is introduced by numerical reasons. 

c) In many cases the effective interaction V^g a ,g, can formally be written 
as a sum of separable terms. 

^V = E^ ph ^V ( 51 ) 

c 

where Q c are single particle operators characterized by the channel index c. 
As discussed in Appendix |Aj this is particular the case for the effective in- 
teraction of the relativistic point coupling model PC-F1 used in the present 
investigation. Working in the channels given by these operators Q c the nu- 
merical effort can be simplified considerably. 

We insert the effective interaction (51) into the Bethe-Salpeter equation 



(39) and introducing the reduced response function: 



Rcc'iuS) - 22 Qa/3 R ^l3a'/3'(^)Qa'l3'i ( 52 ) 



14 



equation (39) turns into the reduced Bethe Salpeter equation 



R cd (u) = R° cd (to) + J2 R cc M V?Rc"d (u). (53) 



which has the same formal solution as given in Eq. (41). In all cases, where 



one has a continuous channel index c, as for instance the radial coordinate r, 



this is an integral equation. In Eq. (53) the interaction V^ ph is diagonal with 



respect to the channel index c. This is not always the case. However, as we 
shall see in Appendix [A| the relativistic interaction PC-F1 can be expressed 
to a large extent in this way. We have to allow only in specific cases also 
for non-diagonal interactions V^J 1 , as for instance in the case of the Coulomb 
force or in the case of derivative terms. This is a rather simple extension of 
the present method and therefore, for the sake of simplicity, we will restrict 
ourselves in the following to an interaction diagonal in the cannel index c. If 
the external operator F in Eq. d34|) can be expressed by the operators Q c as 



F = J2fcQc (54) 

c 

we finally obtain the strength function as 

S(u) = --lmR FF = --ImVp^H/,. (55) 

7T 7T z — ' 

cc' 

If F cannot be expressed in terms of the operators Q c we obtain Rff from 
the Bethe-Salpeter equation fl39| as 



R FF {") = R°ff(") +J2 R °Fc(")Vc(1 - R\u)V» h yjR%{u). (56) 



4 Treatment of the continuum. 

As we have briefly discussed earlier, a proper treatment of the continuum 
is not possible by using a discrete basis, because one needs a tremendously 
large number of p/i-states to fill up the continuum with. Instead, it can only 
be properly taken into account if one makes use of the more flexible linear 
response formalism in an appropriate channel space. 



Starting from Eq. (52 ) for the reduced response function and using Eq. (47) 
we derive the following expression for the reduced free response, which de- 
pends only on the energy u and the channel indices c, d: 
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E 

ph 



(h\Qi\p)(p\Qc\h) ( P \Qt\h)(h\Q c/ \ P ) 



to - e p + e h 



U + Er 



£h 



(57) 



where h stands for occupied (hole) and p for unoccupied (particle) states. 
It is easy to show that the sum over p can be safely extended to run over 
the full space, since terms of the form ^2 hh , vanish due to the cancellation of 
forward and backward going parts. Using completeness we obtain: 



Rl,(r,r':uj) 



J2(h\Qt iQe-Q c rQi\h) 

^ u + e h -h lo - e h + h 

J2(h\QtG(uj + £ h )Q c , + Q & G(-u + e h )Q+\h).(5&) 



Here, h is the Dirac hamiltonian (|22j) and G(E) = 1/(E — h) is the 
corresponding single particle Green's function. 

In this work we use relativistic zero range forces, thus it is appropriate 
to work in coordinate space. The method described in the following is a 
relativistic generalization of the method introduced by Bertsch et al [31] for 
non-relativistic zero range forces. In this case we solve the response equation 
in r-space, which is considerably simpler than the method introduced in 
Refs. [H] for finite range forces. 

In coordinate representation the indices a, (3,... in Eq. (35) are abbre- 
viations for the "coordinates" 1 = (ri, d\, Si, ti), where s is the spin, t the 
isospin coordinate, and d = 1, 2 labels large and small components. Starting 
from the energy density functional ([7]) we find the effective interaction in 
Eq. (Eol to be of the form (JBTl): 



oo 

V ph (l,2) = J2j dr v c (r) Q ] } 2 \i 



(59) 



with the local channel operators Q c {r) defined by 



Q { c\r) 



8{r-n) (1) 
Id 



r (1) 



(60) 



where we distinguish the "coordinates" abbreviated by the upper index (1) 



and the channel index (r, c) used in Eq. (51). Due to this r-dependance, the 
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dimension of the matrix R^j(r, r'; to) in the numerical applications will be 
the number of r-mesh points times eight, which represents the number of the 
covariant channels c, given in Table [6] of the Appendix [A| This implies that 
all scalar, longitudinal, and transverse modes (isoscalar and isovector) are 



fully included and mixed by the matrix inversion of Eq.(41) 



This channel index has now a continuous part given by the radial co- 
ordinate r and a discrete part characterized by the quantum numbers c = 
(D, S, L,T) where the Dirac index D runs over three 2x2 matrices 7^ = 
70, 1,75 defined in Eq. (80), S = 0, 1 is the spin, L the orbital angular mo- 
mentum and T = 0, 1 the isospin. Further details are given in Appendix [A| 



Inserting the channel operators (60) into Eqs. (52) and (58) we obtain 
the reduced free response function: 

^ c ° c ,(r,r» = ^{g^Wb^^r'^ + e^lMr')) (61) 

hK 

+ Q* h c K QiMr')\l D >G K (r', r; -u + e h )j+\h(r)) 

The sum runs over all the occupied states (hole) states h with the 2-dimensional 
radial Dirac spinor (h(r)\ = (/£(r) g*h{ r )) m Eq. (33) and over all the quan- 
tum numbers k = (Ij) compatible with the selection rules in the reduced 
angular and isospin matrix elements 



Q 



er e WII [ a s c Y Lc ]j\\K), 



(62) 



where e*r c = 1 in the isoscalar channel (T c = 0) and ex — ±1 (for protons or 
neutrons) in the isovector channel (T c =1). The reduced matrix elements of 
the operator [o"5 c ^l c ] j contain integrations over the orientation angles VL and 
sums over the spin indices. The matrix elements of the form {h\^fDG{E) r )D i \h) 
depend on r and r' and are obtained by summing over the Dirac indices 
d = 1, 2 for large and small components. 

The Green's function G K (r,r' , E) describes the propagation of a particle 
with the energy E and the quantum numbers K from r to r'. It can either 
be calculated by spectral or non-spectral methods. In the spectral represen- 
tation [32] it is obtained discrete sum 



G K (r,r';E) = J2 



\n(r)) (n(r')\ 

E — E n 



(63) 
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over a complete set of eigenstates |n(r)) of the radial Dirac equation (33) 
with the quantum number k using box boundary conditions (or an oscilla- 
tor expansion). In this case the continuum is discretized, in correspondence 
to the bound states inside the potential. In principle, the radial quantum 
number n runs over the whole single particle basis characterized by the an- 
gular quantum number k, but one can show that this is identical to summing 
only over the unoccupied states, since the hole-hole pairs in Eq. (62) are not 
contributing, due to the cancellation between forward and backward going 
part. Furthermore, because of the no-sea approximation the states in the 
Dirac empty and therefore the sum over n in Eq. d63J) has also to be 

extended over the negative energy states. This corresponds to the sum over 
the a/i-components discussed in the introduction. In practical applications 
one has to restrict this infinite set by a finite sum introducing an upper limit 
e P — €h < E^ t in energy for the particle states p above the Fermi surface and 
a lower limit e a — > —E^ t for the negative energy solutions a is introduced 
in order to make the - otherwise infinite - sum, tractable. This leads to a 
discretized spectrum. 

In the spectral representation the response function 71° (u) has poles at 
the p/i-energies u = ±(e p — £h) an d the full response function H{lo) has poles 
at the eigenenergies f2 M of the RPA-equation (46) in the same restricted 
space. For real frequencies u it is purely real, and therefore the strength 
function vanishes everywhere apart from these poles. For complex energies 
u + iA/2, however, these poles are shifted from the real axis and one obtains 
a continuous spectrum, with the phenomenological width A. This procedure 
yields identical results as the diagonalization of the RPA-matrix in (46) along 
with a subsequent folding with a Lorentzian as discussed in Eq. (50). 



In the non-spectral or continuum approach [31 J the single particle Green's 
function is constructed at each energy from two linearly independent solu- 
tions of the Schroedinger equation with different boundary conditions at 
r = and at r — > oo. In the relativistic case the Dirac-equation in r-space 
depending on the quantum number k is a two-dimensional equation and 
therefore the corresponding single particle Green's function is a 2x2 matrix. 
Using the bracket notation of Dirac for the 2-dimensional spinors we can 
write 1 331: 



r< : i n J \ w K(r))(uJr ) for r > r , N 

GJr,r;E) = < 7 ,( . , (64) 

1 \u K (r)) (w* K (r')\ for r < r' v ' 
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where u(r) and w(r) are two independent Dirac spinors 

MOHgg), KW>=g3) (65) 

normalized in such a way that the Wronskian 

W = /«,(r)^(r)-^(r)/ u (r), (66) 

which is independent of r, is normalized to unity. The solution u R {r) is 
regular at the origin and the solution w K (r) fulfills outgoing wave boundary 
conditions [31] . Further details are given in Appendix |Bj 

Provided that the free response function TZ Q C d (r, r'; lu) has been properly 
derived, we are able to solve the reduced Bethe-Salpeter equation (53) 

Tl C)C >(r,r';uj) = 1Z° c ,(r,r';uj) (67) 

+ V f dr" Tll tC n (r, r"; u) Vc "^ 2 - TZ C »^ (r", r';uj). 

c" J ° T 

where the index c" runs over the various discrete channels given in Table |6j 
Finally the strength function is obtained as: 

S(u) = IhyIZff 

7T 

oo 

\ lm JJ drdr'F*(r)TZ ccl (r, r'; u)F d {r'). (68) 



The sum rules are defined as moments of the strength function S(u): 

POO 

m k = u k S(u) duo. (69) 



They are helpful to characterize the spectral distribution of the oscillator 
strength. In particular they allow us to define the centroid energy by the 
ratio 

E c = (70) 
m 

This quantity can be compared directly with experimental values. Of course, 
in most experiments only a restricted energy range is accessible and therefore 
one also has to restrict the integration in Eq. (69 ) to the same energy window. 
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Other important quantities are transition densities in various channels c 
with respect to the operator F 

/DO 
dr"R cd {ry-uj)F d {r') (71) 

a - 

as for instance the neutron and proton transition densities: 

8p(r)n,p = Sp T =o(r; u) ± Sp T=1 (r; u) (72) 

5 Applications 

In the previous section we briefly described how conventional RPA methods 
treat the continuum part of the spectrum through the introduction of a 
potential "wall" far from the nucleus. In the credit side of this approach, 
general properties of collective excitations can be very well reproduced, either 
by using finite range or point coupling interactions (Niksic et.al. |21j). Since 
CRPA can treat the coupling to the continuum exactly, it is of interest to 
see how well this model does in reproducing the properties of excited state 
in finite nuclei, in particular the giant resonances. 

The most prominent resonances are the Isoscalar Monopole Resonance 
(ISGMR), which is a breathing of the nucleus as a whole, the Isovector 
Dipole Resonance (IVGDR) which corresponds to a collective excitation of 
the proton against the neutron density, and Isoscalar Quadrupole Resonance 
(ISGQR). In addition we have the Isoscalar Dipole Resonance (ISGDR) re- 
vealing the spurious state corresponding to a translational motion of the 
nucleus. These modes show up in an energy range of 10 — 30 MeV and they 
exhaust a major portion of the corresponding sum rules. In the next sections 
we investigate the ISGMR, the IVGDR and the ISGDR in more detail. 

Numerical details 

In the following, we perform several calculations using the relativistic con- 
tinuum RPA approach in r-space with Point Coupling forces [27] . We select 
the doubly magic nuclei 16 0, 40 Ca, 132 Sn and 208 Pb to investigate how the 
collective excitation phenomena depend on an exact coupling to the contin- 
uum. 
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In a first step, the ground state of the nucleus is determined by solving 
the self-consistent RMF equations (33) for the parameter set PC-F1 given in 
Table [1} The method we are using is a fourth order Runge-Kutta in r-space 
(Dirac-mesh) where nucleons move in a spherical box with radius Ro = 15 
fm and with a mesh size dn = 0.05 fm. 

Using the single particle wave functions and the corresponding energies 
of this static solution, we determine the free response 1Z° of Eq. (62) in the 
same box radius but using a wider mesh in r-space (response-mesh). The size 
cIr of this mesh depends on the excitation mode; for the monopole modes 
we use dn = 0.15 fm, while for the dipole a larger interval da = 0.30 fm is 
sufficient. Then we solve the Bethe salpeter equation (68) to get the strength 
distribution S(u). 

At the same time, we perform similar calculations using the discrete RPA 
approach, where the continuum is not treated exactly, aiming of course to a 
more precise comparison with the CRPA results. For those calculations, an 
energy cut-off is necessary, so that a feasible diagonalization is achieved. In 
particular, we have used an energy cut-off \e p - e h \ < E^ut = 300 MeV for the 
configurations with particles above the Fermi sea and \e a — th\ < E^ t = 1500 
MeV for configurations with anti-particles in the Dirac sea. 



Isoscalar Giant Monopole Resonances 

Results for the isoscalar monopole strength distribution are attainable, once 
the corresponding external field 

*£b° = X>? (73) 

i 

is used. In this case, the classical energy weighted sum rule mi (E0) becomes: 

1 ft 2 2ft 2 

mi(E0) = -([F,[T,F]}) = ^(V 2 F) = — (r 2 ). (74) 

The doubly magic spherical nucleus 208 Pb is a particularly good example 
in perform our calculations, since it has been used in the literature to test 
numerous nuclear structure models in the past, in particular applications of 
the random phase approximation [351 E2J [23], [36] . 

In Fig. [I] we show the ISGMR strength distribution obtained by contin- 
uum RPA (full red line) and compare it with the discrete B(E0) values (blue) 
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Figure 1: (Color online) (a) The isoscalar monopole spectrum in 208 Pb, cal- 
culated with the parameter set PC-F1. The red curve corresponds to the 
strength distribution (units on the l.h.s.) obtained by a non-spectral repre- 
sentation without smearing (A = 0), the blue lines give the discrete B(E0)- 
values (units on the r.h.s.) obtained by the spectral representation with the 
same force. The black arrow indicates the experimental centroid energy of 
the resonance [37] . (b) the neutron and proton transition densities at the 
peak with the energy E = 14.40 MeV. 



obtained by the spectral representation of the response function for the same 
parameter set PC-F1 [27j . 

Using the CRPA approach, we find for the calculated centroid energy 
defined in Eq. (70) that rrii/mo = 14.40 MeV, which is rather close to the 
result mi/rriQ = 14.17 MeV deduced from discrete RPA calculations as well 
as to the experimental value mi/m = 13.96 ± 0.2 MeV [37J. 

In those two methods, no additional smearing A = has been used. This 
means that the observed width of the continuum RPA strength corresponds 
entirely to the escape width which in the Pb region is very small, due to the 
relatively high Coulomb and centrifugal barriers in this heavy nucleus. In 
contrast, discrete RPA provides no width at all. Otherwise, the agreement 
of these two methods in this nucleus is excellent. 
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Figure 2: (Color online) The isoscalar monopole strength distribution for 
doubly magic nuclei (a) in 16 0, (b) in 40 Ca, and (c) in 132 Sn. Details are the 
same as in the panel (a) of Fig. [T] 



In the panel (b) of Fig. [TJ we give the neutron and proton transition 
densities at the peak energy, as it is calculated in Eq. (72). They emphasize 



the collective character of the isoscalar breathing mode extended over the 
entire interior of the nucleus with neutrons and protons always in phase. 

In addition, the energy weighted sum rule obtained in CRPA using Eq. 
(69) is mi(EO) = 5.448 ■ 10 5 [MeV-fm 4 ]. This result is in excellent agreement 
with the DRPA calculation mi (SO) = 5.446 ■ 10 5 [MeV-fm 4 ] as well as the 
classical value m^EO) = AAh/2m(r 2 ) = 5.453 ■ 10 5 [MeV-fm 4 ]. This shows 
that the results obtained in the literature by relativistic RPA calculations 
using the spectral method are very reliable for such heavy nuclei [TTl IE]- 

In Fig. [2] we show the EO strength distributions for the lighter doubly 
magic nuclei 16 0, 40 Ca, and 132 Sn. As in Fig. [TJ the smearing parameter A is 
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Figure 3: (Color online) (a) The ISGMR centroid energies as a function of 
the mass number, (b) The experimental and theoretical width of the ISGMR 
as a function of the mass number. Details are given in the text. 

zero, but now the escape width is considerably larger for these nuclei. Fig. [3] 
summarizes the results for the isoscalar monopole strength distributions as 
a function of the mass number A. In panel (a), we plot the centroid energies 
of both continuum RPA (red dots) and discrete RPA (blue dots), together 
with the experimental centroid energies taken from Ref. |37J. We also show 
the phenomenological A-dependence E\- ~ 31.2 A" 1 / 3 + 20.6 A^ 1 ^ by the 
dashed line. It becomes clear that CRPA can successfully reproduce collective 
excitations over the known range of nuclei. 

In panel (b) of Fig. [3] we show the escape width of EO resonances. 
The red values correspond to the full width half maximum (FWHM) of the 
peak, using continuum RPA , while the experimental values are indicated in 
black. The evident disagreement is not surprising, if we consider that only 
lpl/i-configurations are taken into account, i.e. the major part of the width 
resulting from the coupling to more complicated configurations such as 2p2h 
etc. is not described well in this simple RPA approach. It has been shown 
in recent investigations of the coupling to complex configurations within the 
framework of the relativistic time-blocking approximation (RTBA) or the rel- 
ativistic quasiparticle-time-blocking approximation (RQTBA) [21] that such 
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couplings can be taken into account successfully in a fully consistent way 
starting from one density functional E[p\. So far, relativistic investigations 
of this type have been carried out with discrete methods. At present, inves- 
tigations in this direction including the continuum properly go beyond the 
scope of this paper. 



Isovector Giant Dipole Resonances 

Isovector Giant Dipole resonance is the most well studied collective excitation 
and the first to be observed experimentally [SB]- An external electromagnetic 
field of the form: 

N Z Z - 

F L=l r P Y M^ P ) - J E r nY 1M ^n) (75) 

p=l n=l 

causes protons and neutrons to oscillate in opposite phases to each other and 
this leads to a pronounced peak in the photoabsorption cross section. This 
mode has been well studied in many nuclei [39] . 

With the increasing number of experiments in systems far from stability 
and systems with large neutron excess, one has been able to observe also low- 
lying El strength in the area of the neutron emission threshold. It is called 
Pygmy Dipole Resonance PDR and can be interpreted as a collective mode 
with dipole character where the neutron skin oscillates against an isospin 
saturated proton-neutron core. This mode has first been predicted in phe- 
nomenological models [10] exhausting several percent of the electric dipole 
sum rule. In recent years, it has been intensively investigated both on the ex- 
perimental side by the Darmstadt group [2J H2] as well as on the theoretical 
side, using discrete relativistic RPA calculations based on NL3 [43j. 

In Fig. [4] we show in panel (a) the results of the isovector dipole strength 
El in the nucleus 208 Pb using the CRPA approach. The centroid energy at 
13.32 MeV is in excellent agreement with the experimental excitation energy 



E = 13.3 MeV |44| . The energy weighted sum rule (69) is found as mi(El) = 
916.28 [MeVfm 2 ]. This result is in agreement with the DRPA calculation, 
where we obtain mi(El) = 943.32 [MeV-fm 2 ] and as usual somewhat (23.8 
%) larger than the classical Thomas- Reiche-Kuhn sum rule 

9 h 2 NZ 

m TR K = ~r^ ~r = 740.13 [MeV • fin 2 ]. (76) 

47r 2m A 
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Figure 4: (Color online) (a) The isovector dipole strength distribution in 
208 Pb. Details are essentially the same as in the panel (a) of Fig. [TJ However, 
in order to distinguish the continuum (red curve) and the discrete (blue lines) 
calculations we have used here a small smearing parameter A = 10 keV in 
the continuum calculation. The black arrow indicates the theoretical neutron 
emission threshold, (b) transition densities for neutrons and (c) for protons 
at the energy of the PDR (left) and at the GDR (right). 



In addition to the giant dipole resonance a smaller peak appears at the en- 
ergy region of the neutron emission threshold around E ~ 7.5 MeV, that 
corresponds to the pygmy resonance. 

In panel (b) of Fig. [4] we give the transition densities associated the low- 
lying peak at E = 7.66 MeV and the GDR peak at E = 12.9 MeV. The higher 
peak has clearly an isovector character, since the neutrons are oscillating 
against the protons over a large radial range centered at the surface. The 
lower peak shows an isoscalar core, where neutrons and protons oscillate in 
phase and a pure neutron skin moving against the T = core. This is the 
typical behavior of the pygmy mode. 

Closer investigation of pygmy resonances have shown that this mode is 
in the neighborhood of the neutron separation threshold, slightly below for 
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Figure 5: (Color online) The El pygmy resonance (PDR) in the nucleus 
208 Pb. The black arrow indicates the theoretical neutron emission threshold 
at 7.58 MeV. The red dashed lines are obtained by CRPA calculations below 
the threshold. 



small and slightly above for large neutron excess (see for instance Ref. [2"U]). 
It is therefore of particular importance to study this mode with a proper 
treatment of the continuum, since in most of the previous investigations this 
has not been possible. We show in Fig. [5] the details of the PDR in the 
nucleus 208 Pb. Above the theoretical neutron separation threshold which is 





CRPA 


DRPA 


No. 


E 


B(E1) 


E 


B(E1) 


1 


6.90 


0.19 


7.12 


0.23 


2 


7.44 


1.45 


7.46 


2.82 


3 


7.66 


1.11 


7.69 


0.40 


£ 




2.75 




3.45 



Table 2: Energies and B(E1) values for the three most dominant peaks in the 
PDR area around the neutron threshold for the nucleus 208 Pb for continuum 
(CRPA) and discrete (DRPA) calculations. 
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Figure 6: (Color online) The isovector dipole strength distribution in 132 Sn. 
Details are the same as in panel (a) of Fig. |4j 



found at E t ^ = 7.58 MeV (black arrow) we have a continuous red curve 
showing the El strength distribution calculated with CRPA (units at the 
l.h.s) and also few full blue vertical lines that correspond to the discrete 
poles of the DRPA equations (J63j) (units at the r.h.s.) and with length equal 
to the corresponding B(E1) values. 

In the same figure and below the threshold we have in both cases discrete 
lines. The solid blue ones are again the eigen-solutions of the DRPA-equation 
(46 ). The solutions of the CRPA equations lead in this region also to discrete 
poles. We show them by dashed red lines at the pole of the full response 
function. Numerically, the only way to determine the B(E1) values of these 
poles in CRPA is by using very small imaginary parts A — > in the frequency 
uj + il, A and then determining the B(E1) values by simple integration over 
a small interval around this pole. 

By doing that, we finally observe that there are differences in the details 
between the continuum and the discrete RPA calculations close to the neu- 
tron separation threshold. In Table [2] we show for both calculations the three 
most dominant peaks in the area of the PDR around 7.5 MeV. In the discrete 
calculations (DRPA) the strength is concentrated in one peak at E = 7.46 
MeV, whereas in the continuum calculations (CRPA) most of the strength 
in this region is distributed over two peaks, one below the neutron thresh- 
old at E — 7.44 MeV and a sharp resonance slightly above the threshold at 
E = 7.66 MeV. The energy weighted strength in this area is 17.09 [e 2 fm 2 ] 
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Figure 7: (Color online) The El pygmy resonance (PDR) in the nucleus 
132 Sn. Details are the same as in Fig. [Hj The arrow indicates the theoretical 
neutron emission threshold at 22 t h = 7.13 MeV. 



(i.e. 1.86 % of the total sum rule) for CRPA and 26.95 [e 2 fm 2 ] (i.e. 2.85 % 
of the total sum rule) for DRPA. 

In Fig. [6] we show the distribution of the isovector dipole strength in the 
doubly magic nucleus 132 Sn. Again, results using continuum RPA equations 
(red curve) are compared with the solutions obtained from the spectral rep- 
resentation (blue lines). As one can see, there is excellent agreement between 
the two methods, as far as the resonance position and the overall distribution 
is concerned. Moreover, the energy weighted sum rule obtained in CRPA is 
given by mi (221) = 563.60 [MeV-fm 2 ], which is in very good agreement with 
the DRPA calculation mi (El) = 591.02 [MeV-fm 4 ] and 22,9 % larger than 



the Thomas- Reiche-Kuhn sum rule in Eq. ( 76 ) 



In addition, we find that the escape width in this nucleus is considerably 
smaller in the El channel as compared to the E0 channel in Fig. |2j This 
has the following explanation: The selection rules for p/i-excitations with E0 
character is Aj = and no change in parity. It turns out that most of the 
p/i-excitations contributing to the strong peak in the resonance region have 
rather small i values for the particle configurations and therefore a very low 
or no centrifugal barrier. This is different for the El resonance, where one 
has a change in parity and in addition changes of Aj = 0, ±1. In such a case, 
a large part of the contributing p/i-pairs have particles with larger ^-values 
i.e. a strong centrifugal barrier and hence the width becomes smaller. 
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CRPA 


DRPA 


No. 


E 


B(E1) 


E 


B(E1) 


1 


8.11 


0.03 


8.067 


0.037 


2 


8.48 


0.02 


8.186 


1.601 


3 


8.82 


1.44 


8.511 


0.260 


S 




1.490 




1.898 



Table 3: Energies and B(E1) values for the three most dominant peaks in 
the PDR area above the neutron threshold for the nucleus 132 Sn for contin- 
uum (CRPA) and discrete (DRPA) calculations. The units are MeV for the 
energies and [e 2 fm 2 ] for the B(E1) values. More details are given in the text 




E [MeV] 

Figure 8: (Color online) The isovector dipole strength distribution in the 
nuclei 16 (a) and 40 Ca (b). Details are the same as in panel (a) of Fig. 
|4| The theoretical neutron separation energies, indicated by black arrow are 
E thr = 11.33 MeV for 16 and E fhr = 8.91 MeV for 40 Ca. 
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CRPA 


DRPA 


Exp. 


16Q 


20.6279 


21.623 


23.35±0.12 US] 


40 Ca 


18.367 


19.32 


21.76±0.11 [16] 


132 Sn 


14.503 


14.78 




208pj^ 


13.32 


13.23 


13.3±0.10 |H] 



Table 4: Isovector dipole (IVGDR) excitation energies in [MeV] for several 
spherical nuclei, calculated with both continuum and discrete relativistic 
RPA based on the point coupling force PC-F1. 

In Fig. [7] we show the region of the PDR in the doubly magic nucleus 132 Sn. 
As already found in Ref. [20], the theoretical neutron emission threshold at 
E — 7.13 MeV lies much below the area of interest. As before, we calculate 
the B(E1) values of the prominent peaks, for both discrete and continuum 
calculations with the total strength to be in good agreement. In Table [3] 
we show in what extent each level contributes to the total pygmy collective 
state. Finally, the energy weighted strength mi in this area is 13.24 [e 2 fm 2 ] 
(i.e. 2.35 % of the total sum rule) for CRPA and 20.45 [e 2 fm 2 ] (i.e. 3.46 % 
of the total sum rule) for DRPA. 

In Fig.[8]we show the electric dipole strength distribution of the lighter nu- 
clei 16 and 40 Ca. The strength obtained in CRPA calculations (red curves) 
are compared with the B(El)-values resulting from discrete DRPA calcula- 
tions (blue lines). The position of the corresponding peaks and poles with 
large strength are in rather good agreement, as explained in Table |4| We 
find, however, that in the continuum calculations a much larger escape width 
emerges, in particular for the nucleus 16 O. 

Isoscalar Giant Dipole Resonances 

Besides the distribution of the isovector dipole strength which is dominated 
by the IVGDR in many experimental spectra, in recent years there has also 
been considerable interest in measuring the isoscalar dipole strength distri- 
bution [17. '48]. In a similar way, one expects to find the ISGDR, which 
corresponds to a compression wave going through the nucleus along a def- 
inite direction and to learn from such experiments more about the nuclear 



31 



incompressibility. Relativistic calculations based on discrete RPA [T7J [23] 
have shown that the resonance energy of this mode is indeed closely con- 
nected to the incompressibility of nuclear matter. 

Along with this ISGDR resonance built on 3 ^-excitations above 20 MeV, 
calculations based on both relativistic [T7j and non-relativistic RPA ap- 
proaches have revealed a low-lying isoscalar dipole strength in the region 
below and around 10 MeV. Experimental investigations with inelastic scat- 
tering of a-particles at small angles [501 SB] have also found isoscalar dipole 
strength in this region. This strength has been attributed in Ref. [19] to 
an exotic mode of a toroidal motion predicted already in early theoretical 
investigations on multipole expansions of systems with currents [51 J and in- 
vestigated also by semiclassical methods [52] 

On the theoretical point of view, there is further interest in the isoscalar 
dipole mode, characterized by the quantum numbers (J n = 1~,T = 0), be- 
cause it contains the Goldstone mode connected with the violation of trans- 
lational symmetry in the mean field solutions. This mode corresponds to the 
center of mass motion of the entire nucleus. Because of the missing restoring 
force, this mode has vanishing excitation energy. It is one of the essential ad- 
vantages of the RPA approximation, that it preserves translational symmetry 
and therefore it has an eigenvalue at zero energy with the eigenfunction given 
by the p/i-matrix elements of the linear momentum operator. 

Since the ISGDR is expected to be a 3fta>-excitation it is usually associated 
with the external field derived in Ref. [M] 



F l=i = E( r *-^) y ^)> ( 77 ) 

i 

where the factor 77 = |(r 2 ) is used to extract the spurious center of mass 
motion. 

In the upper part of Fig. [9] we display the distribution of the isoscalar 



dipole strength in 208 Pb, calculated with the operator (77) for 77 = 0, that is, 
we take no action for the spurious state. We therefore observe a huge peak 
close to zero energy, which dominates the spectrum and corresponds to the 
spurious translational mode. 

It turns out that the position of this spurious state is an extremely sensi- 
tive object which strongly depends on the numerics of the model. Of course 
the optimal would be to calculate the spurious state at exactly zero energy. 
Therefore this excitation mode presents an ideal benchmark for numerical 
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Figure 9: (Color online) (a) Spurious El isovector strength distribution in 
208 Pb obtained by CRPA calculations with two different values of the radial 
mesh size Sr. (b) the position of the spurious El-state as a function of the 
radial mesh size 



efficiency of the RPA or the linear response equations. Detailed studies have 
shown that the exact separation of the spurious state requires a fully self- 
consistent solution [22]; a fact which was not given in most of the older 
applications with Skyrme or Gogny forces. In many cases, only few of the 
different terms in the residual interaction had been taken into account in 
RPA calculations. 

In addition, the configuration space must be full. Indeed, the discussed 
drawback of the conventional spectral representation in a truncated ph- 
configuration space affects the position of the spurious state. Therefore, 
the convergence to zero eigenvalue of the spurious translational mode occurs 
very slowly and only in extremely large configuration space. In relativistic 
applications this is translated to including also large spectrum in the Dirac 
sea [IS1IH]- As a consequence, in the spectral representation, one has to take 
into account many configuration with particles in the Dirac and holes in the 
Fermi sea, which complicates the numerical applications considerably and 
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inevitably decreases the efficiency of the method. 

Fortunately, using the continuum RPA approach, one is free from such 
constraints and limitations, since the entire configuration space is automat- 
ically included. The results in Fig. [9] obtained with the operator ( 77 ) for 
1] = show clearly the spurious state dominating the entire spectrum (see 
the scale). Its position is not precisely at zero energy, rather it depends on 
the mesh size used for the solution of the continuum response equation (the 
response mesh). In panel (a) of Fig. [9] we present two calculations with dif- 
ferent mesh-sizes, where in panel (b) we show how the spurious state moves 
to zero energy as we use a finer radial interval. For the ideal case of an 
infinitesimal mesh, the strength connected with the spurious state would be 
completely separated from the rest of the spectrum. 

In Fig. 10 we show results obtained with the full operator (77), i.e. with 
rj = |(r 2 ), in a scale increased by three orders of magnitude. Obviously this 
procedure removes the spurious state with high precision. We also did not 
observe any influence of the isoscalar mode in the isovector channel due to 
isospin mixing. In this context we have to remember, that the isospin mixing 
introduced on the mean field level is corrected on the RPA level to a large 
extend 



The main part of the remaining isoscalar dipole spectrum in Fig. 10 



is 
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Low [Me V] 


High [Me V] 


CRPA 


10.97 


25.05 


Hamamoto et al [55] 


~14 


23.4 


Colo et al [56] 


10.9 


23.9 


Vretenar et al. [17J 


10.4 


26. 


Piekarewicz [23] 


~ 8 


24.4 


Shlomo, Sanzhur [57] 


~15 


~25 


Uchida et al. gS] 


12.7 ± 0.2 


22.4 ± 0.5 



Table 5: Self-consistent (relativistic and non-relativistic) RPA calculations 
performed for the ISGDR in 208 Pb, compared with the most recent experi- 
mental data. The two columns refer to the centroid energies of both the low- 
and high-energy sides of the ISGDR mode. 



located at E ~ 25 MeV. This "exotic" mode is best described as a "hydro- 
dynamical density oscillation", in which the volume of the nucleus remains 
constant and the state can be visualized as a compression wave oscillating 
back and forth through the nucleus 



Moreover, Fig. 10 shows an additional mode in the region of 10 — 15 
MeV that exhausts roughly 20% of the total sum rule. This peak does 
not correspond to a compression mode, but as discussed in Ref. [19] rather 
to a kind of toroidal motion. The toroidal dipole mode is understood as a 
transverse zero-sound wave and its experimental observation would invalidate 
the hydrodynamical picture of the nuclear medium, since there is no restoring 
force for such modes in an ideal fluid. 

In conclusion, continuum RPA calculations manage not only to predict 
the existence of the toroidal and the compression mode, but also to achieve a 
reasonable agreement of the corresponding centroid energies to other models 
focusing on the same problem, as well as to recent experimental data [1HI H7J ■ 
In Table |5j these results are presented for the case of the well studied nucleus 

208 pk 
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6 Conclusions 



Starting from a point coupling Lagrangian, we have used the non-spectral 
relativistic RPA approach to examine the corresponding excitation spectra 
and we have compared the results with spectral calculations based on the 
same Lagrangian. This non-spectral method has several advantages. The 
coupling to the continuum is treated consistently using the relativistic single 
particle Green's function at the appropriate energy. In this way, complicated 
sums over unoccupied states are avoided. This is particularly important 
for relativistic applications since the Dirac sea is now automatically treated 
properly and the unphysical transitions from holes in the Fermi sea to par- 
ticles in the Dirac sea is avoided as long as we restrict our investigations to 
positive energies. 

The ground state phenomena are calculated using the same Lagrangian 
by a self-consistent solution of the relativistic mean field equations in r-space. 
The residual particle-hole interaction used in the RPA calculations is derived 
in a fully self-consistent way from the second derivative of the corresponding 
energy density functional. In this way no additional parameters are required 
and one is able to reproduce the collective properties, namely the multipole 
giant resonances for various doubly closed shell spherical nuclei over the entire 
periodic table. 

The calculations are carried out by using a new relativistic continuum 
RPA program for point-coupling models, that includes all the terms in the 
Lagrangian, in particular the two-body interactions with zero range, the den- 
sity dependent parts with all the rearrangement terms, the derivative terms, 
the various current-current terms and the Coulomb interaction. As appli- 
cations the nuclei 16 0, 40 Ca, 90 Zr, 132 Sn, and 208 Pb have been investigated 
proving that a hight level of accuracy is achieved, as compared to the discrete 
methods. Comparing calculations with spectral and non-spectral represen- 
tations of the response function for the same Lagrangian we find, that in 
general the spectra are well reproduced within the spectral approximation, 
if an appropriate phenomenological smearing parameter is used and if a suf- 
ficiently large number of p/i-configurations is taken into account in the latter 
case. We find, however, differences in neighborhood of the neutron thresh- 
old, where the coupling to the continuum is not properly reproduced in the 
spectral method. 

As compared to the discrete case the non-spectral representation has the 
advantage of (i) a precise treatment of the coupling to the continuum and a 
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fully consistent determination of the escape width without a phenomenolog- 
ical smearing parameter, (ii) a faster evaluation of the cross section, because 
one needs for fixed energy only two scattering solutions instead of the thou- 
sands of p/i-configurations in the discrete case and (iii) a proper treatment 
of the Dirac sea without any further a/i-configurations. 

Relativistic CRPA describes very well the position of resonances in doubly 
magic spherical nuclei. Provided that proper pairing correlations are taken 
into account, a similar method can also be applied in open-shell nuclei. This 
requires the development of the relativistic continuum quasiparticle random 
phase approximation (CQRPA). This approach accounts on equal footing for 
the influence of the residual particle-hole (ph) as well as the particle-particle 
(pp) correlations. In analogy to non-relativistic calculations [5E1 EH1 [601 [61] 
this can be achieved on the basis of relativistic CRPA theory developed in this 
manuscript either by treating the pairing correlations in the BCS approach 
for nuclei far from the drip lines where no level in the continuum is occupied, 
or in the Hartree-Bogoliubov approximation valid for all nuclei up to the drip 
line. Investigations in this direction are in progress. 

Of course, the present approach is based on the RPA and includes only 
lpl/i-configurations. Therefore only the escape width of the resonances can 
be reproduced properly. For heavy nuclei the decay width resulting from a 
coupling to more complex configurations is very important. In fact, such cou- 
plings have been introduced successfully in the relativistic scheme using the 
spectral representation in Refs. [24] . On the non-relativistic side, such tech- 
niques have also been used in the context of the non-spectral representation 
without [621 E3] and with [M] pairing. So far, however, fully self-consistent 
relativistic applications including complex configurations with a proper treat- 
ment of the continuum are still missing. 

Helpful discussions with G. Lalazissis, E. Litvinova, T. Niksic, N. Paar, 
V. Tselyaev, and D. Vretenar are gratefully acknowledged. This research has 
been supported the Gesellschaft fiir Schwerionenforschung (GSI), Darmstadt, 
the Bundesministerium fiir Bildung und Forschung, Germany under project 
06 MT 246 and by the DFG cluster of excellence "Origin and Structure of 
the Universe" (www.universe-cluster.de). 
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A The effective interaction in density depen- 
dent point-coupling models 



In Eq. (40) the effective interaction for RPA calculations is defined as the 



second derivative of the energy functional with respect to the density matrix: 

h _ 5 2 E[p] 

V a/3a'/3' ~ t « • \'°) 

In coordinate representation the indices a,f3, . . . are an abbreviation for the 
"coordinates" 1 = (ri,si,di,ti), where s is the spin and t the isospin co- 
ordinate, and d — 1, 2 is the Dirac-index for large and small components. 
Starting from the energy density functional ([7| and neglecting for the mo- 
ment the Coulomb force, we find the density dependent zero range force 

V*\l,2) = ^rWtfCn-raKCrOrJW (79) 

c 

where the vertices T c are 8x8 matrices acting on the indices s, d, t and reflect 
the different covariant structures of the fields including spin and isospin de- 
grees of freedom. We express the 4x4 Dirac matrices as a direct product of 
spin matrices a and 2x2 matrices 70 acting on large and small components 

and the spin matrices <Js=o — 1 an d <Js=i = a n with the spherical coordinates 
of the Pauli spin matrices. In this way we obtain the vertices T c = 7^ x o" S x 
as direct products of 2-dimensional Dirac-, spin- and isospin matrices (see 
also the second column of Table [6]) . 



Finally, in Eq. (79) the quantities v c (r) describe the strengths of all the 
various parts of the interaction derived in a consistent way from the La- 
grangian. The ones derived from the four-fermion terms ^ are constants. 
Furthermore, due to a density dependence of the higher order terms Q as 
well as the corresponding rearrangement terms, v c (r) depends on the static 
density and therefore on the coordinate r. In addition, because of the deriva- 
tive terms ([5]), they also contain Laplace operators. Summarizing, we have: 

c v c (r) = 

scalar: a s + 2p s ps(r) + 37 5 p|(r) + 5 S A 

time-like vector: ay + 3"fvPv( r ) + 

space-like vector: — ay — 7vPy( r ) — 
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In the isovector case the constants as, ay, $s an d $v are replaced by a-rs, 
arv, 8ts arid 5tv- As we see in Table [T] the corresponding values (3ts = 
Its = Itv vanish. 

For spherical nuclei, the densities and currents in the Lagrangian depend 
only on the radial coordinate r. Therefore we expand the 5-function in Eq. 



( 79 ) in terms of spherical harmonics 



S(ri - r 2 ) 



HTi ~ r 2 ) 
r x r 2 



(82) 



Combining spin (S) and orbital (L) degrees of freedom we find by re-coupling 
to total angular momentum J 



1 o 



(i) 



af){Y L {l).Y L {2)) 



£[^]S 1} 



[°sY L ]f 



(83) 



Inserting this expression into Eq. (81) we obtain for the interaction a sum 



(or integral) of separable terms (channels) 



CO 

^(1,2) = E/ dr Qc\r) v c (r) QtW 



(r) 



(84) 



Each channel is characterized by a continuous parameter r and the discrete 
numbers c = (D,S,L,J,T). The corresponding channel operators Q^c\r) 
are local single particle operators 



Q 



(i)/ 



5{r 



a 



(i) 



r (1) 



(85) 



and the upper indices (1) and (2) in Eq. (84) indicate that these operators 
act on the "coordinates" 1 = {r\Q<iS\d\b\) and 2 = (r 2 Q 2 s 2 d 2 t 2 ) . 

The total angular momentum is a good quantum number and for fixed 
J the sum over c in Eq. 84 runs only over specific numbers c = (D, S, L, T) 



determined by the selection rules. We concentrate in this manuscript on 
states with natural parity, i.e. n = (— ) L = (— ) J . Considering that S = for 
the scalar and the time-like vector and that S = 1 for the space-like vector 
we therefore have 

J for S = 
J ± 1 for S = 1 
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c 


r c = 7d <g> a s ® T T 


D 


s 


L 


T 


1 


7o® 1 ® 1 


S 





J 





2 


1® 1 ® 1 


V 





J 





3 


75® a ® 1 


V 


1 


J -I 





4 


7 5 ® cr ® 1 


V 


1 


J + l 





5 


7o® 1 <8> r 3 


s 





J 


1 


6 


1 ® 1 ® r 3 


V 





J 


1 


7 


75® a <g> r 3 


V 


1 


J -I 


1 


8 


75® £7 ® r 3 


V 


1 


J + l 


1 



Table 6: Vertices and quantum numbers of the different channels in Eq. (79) 



Finally we have eight discrete channels. Their quantum numbers are shown 
in Table |6l 



An essential feature of the effective interaction (81) is that it contains 
derivative terms in the form of Laplacians A (retardation effects are ne- 
glected). In spherical coordinates, they contain radial derivatives as well as 
angular derivatives. The latter can be expressed by the angular momentum 
operators acting on spherical harmonics Ylm- Therefore we obtain: 



A 



r 2 *d r —~d r 



L(L 



(86) 



Here the radial derivatives d r and d r act on the right and on the left side 
in Eq. (67), i.e. on 1Z®, c {r'r) and on TZ cc n(r,r"). Since the integration is 
discretized r — >• r n = nh the operator d r is represented by a matrix in 
r-space as for instance by the tree-point formula: 



2h 



(^n',n+l — ^n',n-l)- 



(87) 



This means that the term v c (r) in Eq. (68) is no more diagonal in the 
coordinate r and it must be replaced by a matrix v c (r, r'). 

The term which leads to off-diagonal terms in channel space is the Coulomb 
interaction. It brakes isospin symmetry and therefore it will be described by 
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the general form v cc /(r, r'). In particular, we will have 



^c(l, 2) " r 3 )) (1) ^— d(l - t 3 ))( 2 > 

2 |ri-r 2 | 2 

and the r dependance can be written as: 

= y^v c (r,r')Y L (Q) ■ Y L (tf) 



(88) 



(89) 



with 



v c (r, r') 



Ana r; 



2L + 1 r^ +1 ' 



(90) 



and r< and r> are the smaller and the greater of r and r'. This leads to a 
matrix v cc >(r,r') in Eq. (68) as shown in Table [7j 
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/3f 


T 


(XT 






















1 





+ \vc 








-\v c 













- \vc 








-\vc 


/3f 




















T 





-\vc 








+ \vc 





OCT 



















Table 7: The structure of the channel matrix v 
interaction. 



cc>{r, r 



for the Coulomb 



B The continuum representation for the Green's 
function 

In a non-spectral or continuum approach the relativistic single particle Green's 
function G K (r,r'; E) obeys the equation: 

[E - h K {r)) G K (r, r'; E) = 5(r - r'), (91) 
41 



where h K (r) is the radial Dirac-operator of Eq. (33) depending on the quan- 
tum number k = (Ij). This Green's function can be constructed at each 
energy E from two linearly independent solutions 



fu(r] 
9u(r] 



\w(r)) 



fw(r] 
9w(r] 



\u(r)} = 

(u*(r)\ = (f u (r)g u (r)),(w*(r)\ = (f w (r)g w (r)) 
of the Dirac equation with the same energy E 



(92) 
(93) 

(94) 

but with different boundary conditions. The functions u(r) and w(r) are 
normalized in such a way that the Wronskian is equal to: 



E-h K {r)\ \u{r)) = 0, 



E-h K {r)\ \w(r)) = 0, 



W 



fw{r) f u {r) 
9w(r) g u {r) 



f w (r)g u (r) -g w {r)f u {r) = 1. 



(95) 



Of course these scattering solutions depend on the energy E and on the quan- 
tum number k, i.e. we have \u K (r, E)) and \w K (r, E)). The Dirac-equation in 
r-space is a two-dimensional equation and therefore the corresponding single 
particle Green's function is a 2x2 matrix. Using the bracket notation of 
Dirac for the 2-dimensional spinors and following Ref. [33] we can express 
this Green's function as: 

" \w K (r;E))(u* K (r';E)\ 
\u K {r;E)){w* K {r';E)\ 



G K (r,r';E) 



for r > r' 
for r < r' 



(96) 



with 



G K (r\r;E)=G T K (r,r';E) 

The solution u K (r) is regular at the origin, i.e. following Ref. 
for E > V + S in the limit r — > 0: 



u\r) 



ji(kr) 



^iE-V-S) {kr) l 



(2Z+1)!! 
rjE-V 
M fc(2[+l)!! 



(97) 
we have 

(98) 



with k 2 = (E - V - S)(E - V + S + 2m) > and ji(z) is a spherical 
Bessel function of the first kind. The wave function w R {r) represents at large 
distances for E > an outgoing wave, i.e. we have for r — > oo 

rh\ (kr) 

ikr 



w(r) 



\K\E+2m^i '< 
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\k\ E+2m 



ikr 



(99) 



where h\ (z) is the spherical Hankel function of the first kind and for E < 
an exponentially decaying state, i.e. we have for r — > oo 



«,(r)^( J - 1+1 . W i )e- (100) 



where if 2 = (V - 5 - £)(£ - V + S + 2m) > and and #i+i/ 2 (*) 
are modified spherical Bessel functions [65]. For E < the two scattering 
solutions are both real. This absence of any imaginary term will eventually 



give no contribution to the cross section of Eq. (34). We have to keep in mind, 
however, that at energies that correspond to eigen energies of a bound state, 
the solutions u K (r,E) and w K (r,E) coincide up to a factor, which means 
that the Wronskian vanishes at this energy. This corresponds to a pole in 
the response function on the real energy axis. By adding a small imaginary 
part to the energy E — > E + iA we obtain a sharp peak in the strength 
distribution. 



C The free response function in r-space 

The reduced free response 

^0 M = (h\Qt\ P )(p\Qc\h) _ ( P \Qt\h)(h\Q c ,\p) 
^ oo-e p + e h oo + e p -e h 

depends on the energy E and the channel indices c, d . The operators Q c given 



by Eq. (85) are characterized by the channel index c = (r, DSLT). Each 



single particle matrix element of the form {p\Q c \h) in Eq. (101) separates 
into an angular, an isospin and a radial part. 

(p\Qc\h) = (p\r T \h)(K p \\ [a s Y L )j\\K h ){p\j D \h} r . (102) 

Since we consider in this paper only ph-KPA in the same nucleus, the particle 
states have the same isospin as the hole states and thus the isospin matrix 
element (p|rr|/i) is simply a phase ±1. 

Considering that this channel operator has a 5-function in the radial co- 
ordinate, the radial matrix elements {p\^ c \h) r = {p{ r )\lD\h(r)) then depend 
on r. They are found as sums over the large and small components in the 
radial spinors \h(r)) and \p{r)) for fixed values of r. 



43 



The angular matrix elements depend on the quantum numbers k of par- 
ticle and hole states, and, of course, on the channel quantum numbers S and 
L. In particular, we find for S — 0: 



(mYj\\i'f) = 

while for S = 1, it is 

(mwY L ]j\\i'f} = 



i 



-V- 



47T 



3 J 3' 



-3 



1 + 



- V2 



x 47r 
1 L J 
-1 1 



1 L J 




(103) 



3 J 3' 



3 J 3' 
i 1 i 



2 2 

Using for the angular and isospin part the abbreviation 



(104) 



Qph = ( k p\ \ WsY L ]j\\K h }(p\T T \h}, 



(105) 



we obtain for the reduced response function of Eq. (57) in r-space: 

c' (h\lt\v)r{v\lc'\h)r 



n° cc ,(r,r';uj) = J2{Q; C h Q_ 

ph ^ 



ph 



£h 



n *c n c> {h\lAP)r'ip\lt\ h )i 
^hp^hp~ 



uj + e p -e h 



(106) 

As in Eq. (58) we extend the sum over p over the full space and use com- 
pleteness in the radial wave functions: 

KAry-u) = ^{Q^QthH^liG^ry-u + e^Kr')) 

hn 

+ QZQiMr^lcG^r'^-u + e^j+lhir))}. (107) 



Since the angular matrix elements depend only on the quantum numbers 
k the sum over p is here replaced by a sum over the quantum numbers k, 
which is restricted by the selection rules of the reduced matrix elements 
(105). Having the exact form of the Green's function for the static radial 
Dirac equation (33), one can finally construct the non-spectral or continuum 
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reduced response function (62): 



ft°,(rr» = Y,{^QiaUr^ + ^hi h (r';u; + s h ) (108) 

- Qh%QhKlhw( r l " ~ £ hhth( r '; w - e h ) } for r > r' 
= {QKhQthlhu^ w + e h )ii h {r] u + e h ) 

hK 

~ QTkQLiL^uj -e h )li h {r';u -e h )j for r < r' 
where the Dirac matrix elements depend on the coordinate r: 

j c hw (r;E) = (h\ lc \w(E)) r , (109) 

lUr,E) = (h\lc\u{E)) ri (110) 

rt h (r;E) = (u*(E)\ lc \h) r , (111) 

Y wh (r;E) = (w*(E)\ lc \h) r . (112) 



Using Eq. (97) we find 

TZ° dc {r' ) r-uj)=n cc ,{ry ] u) (113) 

It becomes clear now that the undeniable advantage of the non-spectral ap- 
proach as compared to the spectral one, is the fact that the sum over the 
unoccupied states (particle states) is replaced by a sum over the quantum 
number k, which is restricted by the selection rules for the reduced matrix el- 
ements Q c Kh - For each k, one has to determine only the pairs of the scattering 
wave functions \u) and \w) for the forward and backward term. In particular 
the sum over k does not have to be extended over the states in the Dirac 
sea as in the spectral representation (for details see Ref. j8]). Therefore, not 
only the size of the configuration space is significantly reduced, but, more 
notably, the particle-hole as well as the antiparticle-hole basis is taken into 
account fully and without any approximation. 
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